home *** CD-ROM | disk | FTP | other *** search
/ Aminet 2 / Aminet AMIGA CDROM (1994)(Walnut Creek)[Feb 1994][W.O. 44790-1].iso / Aminet / util / misc / Fudgit233.lha / Source / src / lgamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-16  |  8.2 KB  |  304 lines
  1. #ifndef lint
  2. static char sccsid[] = "@(#)lgamma.c    1.2    (ucb.beef)    3/1/89";
  3. #endif    /* !defined(lint) */
  4. /* 
  5.  * This routine calculates the log(GAMMA) function for a positive real
  6.  *   argument x.  Computation is based on an algorithm outlined in
  7.  *   references 1 and 2.  The program uses rational functions that
  8.  *   theoretically approximate log(GAMMA) to at least 18 significant
  9.  *   decimal digits.  The approximation for x > 12 is from reference
  10.  *   3, while approximations for x < 12.0 are similar to those in
  11.  *   reference 1, but are unpublished.  The accuracy achieved depends
  12.  *   on the arithmetic system, the compiler, the intrinsic functions,
  13.  *   and proper selection of the machine-dependent constants.
  14.  * 
  15.  **********************************************************************
  16.  **********************************************************************
  17.  * 
  18.  * Explanation of machine-dependent constants
  19.  * 
  20.  * beta   - radix for the floating-point representation
  21.  * maxexp - the smallest positive power of beta that overflows
  22.  * XBIG   - largest argument for which LN(GAMMA(x)) is representable
  23.  *          in the machine, i.e., the solution to the equation
  24.  *                  LN(GAMMA(XBIG)) = beta**maxexp
  25.  * XINF   - largest machine representable floating-point number;
  26.  *          approximately beta**maxexp.
  27.  * EPS    - The smallest positive floating-point number such that
  28.  *          1.0+EPS > 1.0
  29.  * FRTBIG - Rough estimate of the fourth root of XBIG
  30.  * 
  31.  * 
  32.  *     Approximate values for some important machines are:
  33.  * 
  34.  *                           beta      maxexp         XBIG
  35.  * 
  36.  * CRAY-1        (S.P.)        2        8191       9.62E+2461
  37.  * Cyber 180/855
  38.  *   under NOS   (S.P.)        2        1070       1.72E+319
  39.  * IEEE (IBM/XT,
  40.  *   SUN, etc.)  (S.P.)        2         128       4.08E+36
  41.  * IEEE (IBM/XT,
  42.  *   SUN, etc.)  (D.P.)        2        1024       2.55D+305
  43.  * IBM 3033      (D.P.)       16          63       4.29D+73
  44.  * VAX D-Format  (D.P.)        2         127       2.05D+36
  45.  * VAX G-Format  (D.P.)        2        1023       1.28D+305
  46.  *
  47.  * 
  48.  *                           XINF        EPS        FRTBIG
  49.  * 
  50.  * CRAY-1        (S.P.)   5.45E+2465   7.11E-15    3.13E+615
  51.  * Cyber 180/855
  52.  *   under NOS   (S.P.)   1.26E+322    3.55E-15    6.44E+79
  53.  * IEEE (IBM/XT,
  54.  *   SUN, etc.)  (S.P.)   3.40E+38     1.19E-7     1.42E+9
  55.  * IEEE (IBM/XT,
  56.  *   SUN, etc.)  (D.P.)   1.79D+308    2.22D-16    2.25D+76
  57.  * IBM 3033      (D.P.)   7.23D+75     2.22D-16    2.56D+18
  58.  * VAX D-Format  (D.P.)   1.70D+38     1.39D-17    1.20D+9
  59.  * VAX G-Format  (D.P.)   8.98D+307    1.11D-16    1.89D+76
  60.  * 
  61.  ***************************************************************
  62.  ***************************************************************
  63.  * 
  64.  * Error returns
  65.  * 
  66.  *  The program returns the value XINF for x <= 0.0 or when
  67.  *     overflow would occur.  The computation is believed to 
  68.  *     be free of underflow and overflow.
  69.  * 
  70.  * 
  71.  * Intrinsic functions required are:
  72.  * 
  73.  *      log
  74.  * 
  75.  * 
  76.  * References:
  77.  * 
  78.  *  1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for
  79.  *     the Natural Logarithm of the Gamma Function,' Math. Comp. 21,
  80.  *     1967, pp. 198-203.
  81.  * 
  82.  *  2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May,
  83.  *     1969.
  84.  * 
  85.  *  3) Hart, Et. Al., Computer Approximations, Wiley and sons, New
  86.  *     York, 1968.
  87.  * 
  88.  * 
  89.  *  Authors: W. J. Cody and L. Stoltz
  90.  *           Argonne National Laboratory
  91.  * 
  92.  *  Latest modification: June 16, 1988
  93.  */
  94. #if defined(__STDC__) || defined(__GNUC__)
  95. extern double log(double);
  96. #else
  97. extern double log();
  98. #endif
  99.                     /* Machine dependent parameters */
  100. #if defined(vax) || defined(tahoe)
  101. #define    XBIG    (double)2.05e36
  102. #define    XINF    (double)1.70e38
  103. #define    EPS    (double)1.39e-17
  104. #define    FRTBIG    (double)1.20e9
  105. #else    /* defined(vax) || defined(tahoe) */
  106. #define    XBIG    (double)2.55e305
  107. #define    XINF    (double)1.79e308
  108. #define    EPS    (double)2.22e-16
  109. #define    FRTBIG    (double)2.25e76
  110. #endif    /* defined(vax) || defined(tahoe) */
  111.                     /* Mathematical constants */
  112. #define    ONE    (double)1
  113. #define    HALF    (double)0.5
  114. #define    TWELVE    (double)12
  115. #define    ZERO    (double)0
  116. #define    FOUR    (double)4
  117. #define    THRHAL    (double)1.5
  118. #define    TWO    (double)2
  119. #define    PNT68    (double)0.6796875
  120. #define    SQRTPI    (double)0.9189385332046727417803297
  121.  
  122. /*
  123.  * Numerator and denominator coefficients for rational minimax Approximation
  124.  * over (0.5,1.5).
  125.  */
  126. static double D1 = -5.772156649015328605195174e-1;
  127. static double P1[] = {
  128.     4.945235359296727046734888e0,
  129.     2.018112620856775083915565e2,
  130.     2.290838373831346393026739e3,
  131.     1.131967205903380828685045e4,
  132.     2.855724635671635335736389e4,
  133.     3.848496228443793359990269e4,
  134.     2.637748787624195437963534e4,
  135.     7.225813979700288197698961e3,
  136. };
  137. static double Q1[] = {
  138.     6.748212550303777196073036e1,
  139.     1.113332393857199323513008e3,
  140.     7.738757056935398733233834e3,
  141.     2.763987074403340708898585e4,
  142.     5.499310206226157329794414e4,
  143.     6.161122180066002127833352e4,
  144.     3.635127591501940507276287e4,
  145.     8.785536302431013170870835e3,
  146. };
  147.  
  148. /*
  149.  * Numerator and denominator coefficients for rational minimax Approximation
  150.  * over (1.5,4.0).
  151.  */
  152. static double D2 = 4.227843350984671393993777e-1;
  153. static double P2[] = {
  154.     4.974607845568932035012064e0,
  155.     5.424138599891070494101986e2,
  156.     1.550693864978364947665077e4,
  157.     1.847932904445632425417223e5,
  158.     1.088204769468828767498470e6,
  159.     3.338152967987029735917223e6,
  160.     5.106661678927352456275255e6,
  161.     3.074109054850539556250927e6,
  162. };
  163. static double Q2[] = {
  164.     1.830328399370592604055942e2,
  165.     7.765049321445005871323047e3,
  166.     1.331903827966074194402448e5,
  167.     1.136705821321969608938755e6,
  168.     5.267964117437946917577538e6,
  169.     1.346701454311101692290052e7,
  170.     1.782736530353274213975932e7,
  171.     9.533095591844353613395747e6,
  172. };
  173.  
  174. /*
  175.  * Numerator and denominator coefficients for rational minimax Approximation
  176.  * over (4.0,12.0).
  177.  */
  178. static double D4 = 1.791759469228055000094023e0;
  179. static double P4[] = {
  180.     1.474502166059939948905062e4,
  181.     2.426813369486704502836312e6,
  182.     1.214755574045093227939592e8,
  183.     2.663432449630976949898078e9,
  184.     2.940378956634553899906876e10,
  185.     1.702665737765398868392998e11,
  186.     4.926125793377430887588120e11,
  187.     5.606251856223951465078242e11,
  188. };
  189. static double Q4[] = {
  190.     2.690530175870899333379843e3,
  191.     6.393885654300092398984238e5,
  192.     4.135599930241388052042842e7,
  193.     1.120872109616147941376570e9,
  194.     1.488613728678813811542398e10,
  195.     1.016803586272438228077304e11,
  196.     3.417476345507377132798597e11,
  197.     4.463158187419713286462081e11,
  198. };
  199.  
  200. /*
  201.  * Coefficients for minimax approximation over (12, INF).
  202.  */
  203. static double C[] = {
  204.     -1.910444077728e-03,
  205.      8.4171387781295e-04,
  206.     -5.952379913043012e-04,
  207.      7.93650793500350248e-04,
  208.     -2.777777777777681622553e-03,
  209.      8.333333333333333331554247e-02,
  210.      5.7083835261e-03,
  211. };
  212.  
  213. double
  214. #if defined(__STDC__) || defined(__GNUC__)
  215. lgamma(double x)
  216. #else
  217. lgamma(x)
  218. double x;
  219. #endif
  220. {
  221.     register i;
  222.     double res,corr,xden,xnum,dtmp;
  223.  
  224. #define    y x
  225.     if (y > ZERO && y <= XBIG) {
  226.         if (y <= EPS)
  227.             res = -log(y);
  228.         else if (y <= THRHAL) {            /*  EPS < x <= 1.5 */
  229. #define    xm1 dtmp
  230.             if (y < PNT68) {
  231.                 corr = -log(y);
  232.                 xm1 = y;
  233.             }
  234.             else {
  235.                 corr = ZERO;
  236.                 xm1 = y-HALF; xm1 -= HALF;
  237.             }
  238.             if (y <= HALF || y >= PNT68) {
  239.                 xden = ONE;
  240.                 xnum = ZERO;
  241.                 for (i = 0; i <= 7; i++) {
  242.                     xnum = xnum*xm1+P1[i];
  243.                     xden = xden*xm1+Q1[i];
  244.                 }
  245.                 res = xnum/xden; res = corr+xm1*(D1+xm1*res);
  246. #undef    xm1
  247.             }
  248.             else {
  249. #define    xm2 dtmp
  250.                 xm2 = y-HALF; xm2 -= HALF;
  251.                 xden = ONE;
  252.                 xnum = ZERO;
  253.                 for (i = 0; i <= 7; i++) {
  254.                     xnum = xnum*xm2+P2[i];
  255.                     xden = xden*xm2+Q2[i];
  256.                 }
  257.                 res = xnum/xden; res = corr+xm2*(D2+xm2*res);
  258.             }
  259.         }
  260.         else if (y <= FOUR) {            /*  1.5 < x <= 4.0 */
  261.             xm2 = y-TWO;
  262.             xden = ONE;
  263.             xnum = ZERO;
  264.             for (i = 0; i <= 7; i++) {
  265.                 xnum = xnum*xm2+P2[i];
  266.                 xden = xden*xm2+Q2[i];
  267.             }
  268.             res = xnum/xden; res = xm2*(D2+xm2*res);
  269. #undef    xm2
  270.         }
  271.         else if (y <= TWELVE) {            /*  4.0 < x <= 12.0 */
  272. #define    xm4 dtmp
  273.             xm4 = y-FOUR;
  274.             xden = -ONE;
  275.             xnum = ZERO;
  276.             for (i = 0; i <= 7; i++) {
  277.                 xnum = xnum*xm4+P4[i];
  278.                 xden = xden*xm4+Q4[i];
  279.             }
  280.             res = xnum/xden; res = D4+xm4*res;
  281. #undef    xm4
  282.         }
  283.         else {                    /*  x >= 12.0 */
  284.             res = ZERO;
  285.             if (y <= FRTBIG) {
  286.                 res = C[6];
  287. #define    ysq dtmp
  288.                 ysq = y*y;
  289.                 for (i = 0; i <= 5; i++)
  290.                     res = res/ysq+C[i];
  291. #define    ysq dtmp
  292.             }
  293.             res /= y;
  294.             corr = log(y);
  295.             res += SQRTPI; res -= HALF*corr;
  296.             res += y*(corr-ONE);
  297.         }
  298. #undef    y
  299.     }
  300.     else                    /*  Return for bad arguments */
  301.         res = XINF;
  302.     return res;
  303. }
  304.